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We consider static cosmological solutions along with their stability properties in the framework 
of a recently proposed theory of massive gravity. We show that the modification introduced in the 
cosmological equations leads to several new solutions, only sourced by a perfect fluid, generalizing 
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numerical analysis, we show that the found solutions can be either neutrally stable or unstable 
against spatially homogeneous and isotropic perturbations. 
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I. INTRODUCTION 

The exact solution of Einstein's equations known as 
the Einstein Static (ES) Universe is a static closed 
Friedmann-Robertson- Walker model sourced by a per- 
fect fluid and a cosmological constant (see [1]). Its sta- 
bility properties have been widely investigated. The ES 
Universe is unstable to homogeneous perturbations [2], 
moreover it is always neutrally stable against small inho- 
mogeneous vector and tensor perturbations and neutrally 
stable against adiabatic scalar inhomogeneities with high 
enough sound speed [3]- Furthermore, the ES Universe 
was recently shown to be unstable to Bianchi type-IX 
spatially homogeneous perturbations in the presence of 
nontilted and tilted perfect fluids with p -t- 3P > [5] 
and for several kinds of matter flelds sources (see [1] and 
references therein). 

The renewed interest in the ES Universe, besides its 
historical importance, comes from the Emergent Universe 
scenario [Ij , an inflationary cosmological model in which 
it plays a crucial role as initial state. This model, in 
turn, suffers from a flne-tuning problem which is amelio- 
rated when modifications to the cosmological equations 
of GR are present. For this reason, analogous solutions 
have been considered in several different modified grav- 
ity models [8] and quantum gravity models [9Hl4]. In- 
deed, when dealing with modified cosmological equations, 
many new static solutions are present, whose stability 
properties, depending on the details of the single theory 
or family of theories taken into account, are substantially 
different from those of the classical ES solution of GR. In 
particular, neutrally stable solutions are present thus the 
fine-tuning problem is ameliorated but then a mechanism 
is needed to get out of the phase of infinite expansions 
and recoUapses and to trigger the expanding phase of the 
Universe [7]. 

Here we study static cosmological solutions in the 



framework of a covariant Massive Gravity (MG) model 
recently proposed in |151 116| . In order to construct a con- 
sistent theory, nonlinear terms should be tuned to remove 
order by order the negative energy state in the spectrum 
|17) . The model under investigation follows from a pro- 
cedure originally outlined in [THl [H] and has been found 
not to show ghosts at the complete nonlinear level with 
an arbitrary reference metric (501 HI] ■ 

The considered theory exploits several remarkable fea- 
tures. Indeed the graviton mass typically manifests itself 
on cosmological scales at late times thus providing a nat- 
ural explanation of the presently observed accelerating 
phase [22]. Moreover, the theory allows for exotic so- 
lutions in which the contribution of the graviton mass 
affects the dynamics at early times. Indeed, in contrast 
with GR where, in order to have static solutions, a cos- 
mological constant term and a positive curvature term 
are needed in addition to a suitable perfect fluid source 
term, we find that in the considered MG theory it is pos- 
sible to have static cosmological solutions only sourced 
by a perfect fluid. These solutions can be either unstable 
or neutrally stable and they exist even for spatially flat 
(i.e. /C = 0) cosmological models. 

This paper is structured as follows. In Sec. |TT| the non- 
linear MG framework considered in this work is shortly 
described and the modifled Friedmann equations under 
investigation are introduced. In Sec. |III| static cosmo- 
logical solutions are found, a linearized analysis is per- 
formed and the stability properties are discussed in de- 
tails. In Sec. |IV| the dynamics near the fixed points is 
described using numerical integrations. The phase dia- 
grams of the system are drawn both in the (a, a)— plane 
and {H, p)— plane. In Sec. |vj some conclusions are even- 
tually drawn. 



II. COSMOLOGICAL EQUATIONS 
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We consider the theory introduced in [TSl [H] . In the 
formalism afterwards used in [23], the theory we refer 
to is defined on a four-dimensional pseudo-Riemannian 
manifold {A4,g) and the dynamics is determined by the 
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action 



S = 



8nG 



1 



R + ni^U d*x + Sm, 



where G is the Newton gravitational constant and R is 
the Ricci scalar while Sm describes ordinary matter. The 
potential term, coupled through a mass term m, is de- 
fined by 



C3 



2 1 '^pi'po- 



where e^i,po- is the Levi-Civita tensor density, C3 and C4 
are arbitrary dimensionless real constants and 

7^ being defined by the relation 

with a symmetric tensor field. The quantity nig — 
hm/c is called the graviton mass. 

We consider a Robertson- Walker Universe with three- 
dimensional spatial curvature /C = 0,±1, described by 
the line element 



= dt^-a{tf 



dr^ 



1 -/Cr2 



r'^ide'^ +sin{efdcj)^) 



The first Friedmann equation has been written in [24] 
for generic values of the dimensionless constants C3 and 
C4 with the constraints imposed by Bianchi identities; it 
reads 



TO^ ( 4c3 + C4 - 6 + 3C 



3 - 3c3 



C4 



+2>C 



2 C4 + 2C3 — 1 



a 
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where C is an integrating constant. Matter couples min- 
imally to gravity thus its equation of motion is 



3H{p + p)=0 



(1) 



with H — a/a. From now on we assume a constant equa- 
tion of state parameter w thus p ~ wp. 
Moreover, in the subsequent analysis the parameter space 
is reduced to the subset C3 — — C4 since, as found in pS] . 
this is the simplest choice that presents a successful Vein- 
shtein effect in the weak field limit. 

For later purposes it is useful to rewrite Eq.Q as fol- 
lows: 



m 



A. 



A3 
a2 



(2) 



where k = SttG, a/C — >■ a and 

Ai ~ — 3c4 — 6 
A-2, = 3(3 + 2c4) 

A3 = -3(1 + C4). 

The second Friedmann equation reads 



n 



K. . IC m? f A2 „ A3 



--p{l + w) + — 



(3) 



Making use of the Friedmann constraint Eq.([2| one can 
recast Eq.(|3]) as a second order nonlinear differential 
equation in a and its first and second derivatives 



TO 



(l + w)Ai + 



2 + 3w [A 



l + Zw f A3 



which can be easily recast as a proper two-dimensional 
autonomous dynamical system by introducing the vari- 
ables: 

q — a p — a. 
Thus, the system to be considered is the following: 



q = p 



2q 



Ai{l + w)q + A: 



2 + 3w 



(l + 3w) 



A 



3™2 



TO^ + (/C +p2) 



(4) 



(5) 



The dynamics described by the above equations is glob- 
ally Hamiltonian with respect to the symplectic structure 

oj = q^+^'"dq A dp, 

which is singular in g = 0. Indeed 

ix^^ = —dT-L 

where 



9 3 V q q 



(6) 

and ix is the contraction operator with respect to the 
vector field 



X = P-. 



d 



m I A , -, N , 2 + 3w 
— lAi{w + l)q + A2^^ 



^P» + l)(f„r + (K+/)))|, 

which is singular in g = and p = 0. The Hamilton's 
equations read 



q = 



on 

dp ' 



dn 

dq 
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III. STATIC SOLUTIONS AND THEIR 
STABILITY 

By imposing the condition a = a = p = 0, Eq.Q is 
identicaUy satisfied and the system of Eqs.([I])-([3]) reduces 
to an algebraic system in the unknowns a and p. From 
this we get 



a± = 



P± 



6m(l + w)Ai 



m 



M 

a± 



A. 



(7) 
(8) 



with 

f7 = m2(2 + "iwfAl +12(3/C -mM3)(l +Aw +iw'^)Ai. 

Interesting enough, these solutions may exist not only 
for /C = 1 as in GR, indeed the modified cosmological 
equations of MG allow for static solutions also in both 
the /C = case and the /C = — 1 case. 

The stability analysis of the formerly presented solu- 
tions can be easily performed using standard dynamical 
system techniques. It is easy to check that the solutions 
in Eqs.Q and Q are stationary points of the dynami- 
cal system in Ec^s.Q and ([S]). Their stability is readily 
determined by looking at the eigenvalues of the Jacobian 
matrix evaluated at the stationary points. 

Since the solutions we are interested in are static, once 
the eigenvalues of the linearized system are determined 
in general, we can impose p = 0; they then reduce to a 
particularly simple form 



with 



Ai. 



A3 
3 

+/C(1 + "iw)] . 



± 



V2q^ 



(9) 



{1 + 3w) + Ai{l + w)q' 



(10) 



According to the sign of S evaluated at the stationary 
points, one can have either a pair of real eigenvalues with 
opposite signs or a pair of purely imaginary eigenvalues 
with opposite signs. Thus, the fixed points are either 
neutrally stable [27] (i.e. center) when E < or unstable 
(i.e. a saddle) when E > 0. 

For the sake of simplicity we will first consider in de- 
tails the simplest but interesting case JC = which shares 
all the relevant features of the other more involved ones. 



A. Case 1: K. = 

In the case of spatially flat models the solutions in 
Eq.([7]) reduce to the simple expression 



ai,2 



(2c4 + 3)(2 + 3w) ± ^(2c4 + 3)2 + {3w + 2)2 
6(c4 + 2)(l + w) 



(11) 



in which the explicit dependence on C4 has been restored. 
The corresponding expression for the energy density in 
terms of the model parameters can be obtained replacing 
these solutions in Eq.(|8]). The conditions on the param- 
eters for these solutions to exist can be found by simul- 
taneously imposing a.i > and pi > 0. The results are 
reported in Table |I] 



Sol. 


w 


C4 


Sol.l 


w < -1 


C4 / -2,-1 


-1<TO<-I 


C4/ 2nC4< «+V3(l+4»+3»^) 


<™<-i 


C4 < -2 


™ > -I 


C4 < -2 C4 > -1 


S0I.2 


-l<TO<-f 


^ ^ y3{-l-4»-3»2)+3 



TABLE I: Existence conditions for the static solutions of the 
spatially flat model. 



The stability of the solutions is determined by evalu- 
ating the eigenvalues in Eq.(|9| at these fixed points. As 
already observed, the problem simply consists in evalu- 
ating the sign of the function E in Eq.( 10 ). We find that 



the first solution is always unstable of the saddle type 
while the second is always (neutrally) stable, i.e., it is 
a center. In the latter case the linearized analysis is not 
sufficient to ensure that the solution is actually stable, in- 
deed the hypotheses of the Hartman-Grobman theorem 
are not fulfilled being the point nonhyperbolic. 

To further analyze the second solution we use the Lya- 
punov's second method. Let us consider the function H 
as m Eq.([6]) with /C = and define 



V{q,p)=n{q,p)-n{a2,0). 



(12) 



V{q,p) is positive-definite in a neighborhood U around 
the second solution, i.e. 



Viq,p) 
V{q,p) 



> 



for (g = a2,p = 0) (13) 
V(g,p)eC/\(a2,0), (14) 



thus it is a good Lyapunov function candidate. Its time 
derivative is zero in a neighborhood of the fixed point i.e. 



V{q,p)^Hiq,p)^0 V(g,p)eC/, 



(15) 



thus the second solution is proven to be stable. The 
result is also confirmed by numerically integrating the 



fully nonlinear system (see Sec IV I . 

Figs. |l(a)] and 1(b) | show regions in the (04, w)— plane 



in which Sol.l and Sol. 2 are admitted, respectively; the 
former being unstable while the latter being stable. Such 
regions have a nonvanishing intersection. More precisely. 



Fig. 


1(b) 


1(a) 


For 



contained by the region in Fig. 1(a) For C4 and w in 
this region, both static solutions are admitted; otherwise 
either the unstable solution only or no static solutions 
are admitted. 
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by varying the parameter w while keeping the parameter 
C4 fixed. Further details about this feature are discussed 





(a) 



(b) 



(a) 



(b) 



FIG. 1: Stability conditions for the two static solutions of the 
IC = case in terms of the parameters C4 and w. The stability 
properties of the two solutions in the considered region, —3 < 
C4 < 1 and —2<w<l, are always different. The first 
solution (a) is unstable (U, red); the second solution (b) is 
stable (S, blue). 



FIG. 3: Scale factor value of the static solutions (ai,2) as a 
function of (a) C4 and (b) w. The parameters have been arbi- 
trarily chosen (i.e. « = 1, m = 1, lu = 0.08, C4 = 0.6) in order 
to show a typical bifurcation. The upper (blue) curve and 
the lower (red) curve correspond to the stable and unstable 
solution respectively. 




in the following sections (see Sec IV). 



FIG. 2: This figure is obtained by overlapping Figs. 1(a) and 



1(b) it shows a region (U/S, purple color) of the parameters 
space allowing both solutions to exist, the first being unstable 
while the second being stable. 



It is interesting to note that the phase space of the sys- 
tem exploits relevant changes in its qualitative structure 
according to the parameters values. In particular, the 
system undergoes bifurcations that can be singled out, 
for instance, by looking at the eigenvalues characterizing 
the stability of the fixed points. 

In Fig. [3] an example of bifurcation is depicted. 
Fig. 3(a)| is obtained by varying the parameter C4 while 
keeping the parameter w fixed. The value of the scale fac- 
tor corresponding to each static solution varies until the 
fixed points annihilate exploiting a saddle-center bifurca- 
tion. Fig. |3(b)| shows a similar behavior; it is obtained 



B. 



Case 2: K: = 1 



Spatially closed models can be analyzed following the 
formerly described procedure, that is imposing the posi- 
tivity of a and p and then studying the sign of the func- 
tion E in Eq.(10|. The analytical expression of the re- 



sulting ranges in terms of the parameters w and C4 is 
quite cumbersome and not particularly illuminating, for 
this reason we do not report it explicitly. We just re- 
mark that, as in the spatially flat case (Sec. Ill A), at 



most two static solutions are admitted, one being un- 
stable, the other being neutrally stable, in contrast with 
GR where only one unstable solution of the saddle type 
is admitted. 

In the considered range, — 3 < C4 < 1 and — 2 < ui < 1, 
the existence region of the first solution consists of three 
parts while the existence region of the second solution 
is compact (see Fig. |4]). These regions have two inter- 
sections, which means that for values of the parameter 
within such regions, both solutions are admitted simul- 
taneously, otherwise either only one of the two solutions 
is admitted or no solutions are admitted. This also indi- 
cates that in the spatially closed models the phase space 
has a richer structure and several bifurcations can oc- 
cur. For instance, bifurcations diagrams similar to those 
found in the spatially flat case (depicted in Fig. [s]) can 
be drawn. 



C. Case 3: K. = -1 

As in the case of spatially closed models, the analyti- 
cal expression defining the ranges of existence of the two 
solutions for spatially open models is quite cumbersome. 
As in both the previously analyzed cases, the system in 
Eqs. ([4]) and ([5| admits two static solutions, one being 
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-3 -2 -1 1 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 

C4 C4 



FIG. 4: Stability conditions for the IC — 1 case in terms of 
the parameters C4 and w. The regions of existence of the two 
solutions partially overlap but the character of their stability 
is different, the first solution being unstable (U, red), the 
second solution being stable (S, blue). For parameters in the 
two triangular regions (U/S, in purple color) both solutions 
can exist. 



FIG. 5: Stability conditions for the IC = —1 case in terms of 
the parameters C4 and w. The regions of existence of the two 
solutions partially overlap but the character of their stability 
is always different. The first solution is unstable (U, red) and 
the second solution is stable (S, blue). There is a region of the 
parameter space (U/S, purple color) for which both solutions 
can exist. 



unstable, the other being neutrally stable, in contrast 
with GR where only one static solutions is present [IS]. 

In the considered range, — 3 < C4 < 2 and — 2 < w < 1, 
the existence region of the first solution consists of two 
parts while the existence region of the second solution 
is compact (see Fig. [s]). The existence regions have one 
intersection; for C4 and w in this region, both static solu- 
tions are admitted, otherwise either the unstable solution 
only is admitted or no static solution is admitted. The 
result is very similar to that obtained for open models but 
now the system does not admit solutions for C4 smaller 
than « —2.25. This result also suggests that, crossing 
the boundary of the region, where both solutions are ad- 
mitted, bifurcations can occur. For instance, bifurcations 
diagrams similar to those found in the spatially flat case 
(depicted in Fig. Isl) can be easily drawn. 



IV. NUMERICAL INTEGRATION 

In this section we discuss the properties of the consid- 
ered solutions and the phase space of MG by perform- 
ing numerical integrations of the system in Eqs. (|4| and 
([5]). This procedure allows us to check the results of the 
linearized stability analysis for the nonhyperbolic fixed 
points and provides further interesting physical informa- 
tion. For the seek of simplicity, we just consider the 
spatially flat case. 

Fig.[6]shows a region of the (q,^)— plane for the system 



in Eqs. ([4]) and ([5|. The (red) point on the left represents 
the unstable flxed point, the (blue) point on the right rep- 
resents the nonhyperbolic flxed point. Arrows represent 
the orbits obtained evolving initial conditions. The (red) 
continuous curve is a separatrix, it marks the boundary 
of regions where the dynamical behavior of the system 
is different. Initial conditions, belonging to the region 
enclosed by the separatrix, evolve producing closed or- 
bits which remain close to the nonhyperbolic fixed point 
confirming that it is neutrally stable. 

More physical insights can be gained by rewriting the 
systems in terms of different variables, namely H and p. 
To this aim one first has to consider the original systems 
of Eqs. ([T])-(|3]). The Friedmann constraint in Eq. ([2| can 
be locally solved in order to express the scale factor a in 
terms of H and p, then, by substituting in Eq.(|3|, one 
gets a new equation for H which, together with Eq. ([l]), 
is a two-dimensional autonomous dynamical system. 

Fig. [7] shows the phase space portrait for different suit- 
ably chosen values of the model parameters. The (red) 
continuous line is a separatrix enclosing the stable re- 
gion containing the static solution. The closed orbits 
represent cosmological models characterized by an infi- 
nite sequence of bouncing and recollapsing epochs. The 
pictures Figs. 7(a)| - |7(d)| are obtained by varying the 
value of the equation of state parameter w and keeping 
constant the other model parameter, w plays the role 
of bifurcation parameter indeed, according to its values, 
the distance between the fixed points decreases until the 




5 10 15 



FIG. 6: Phase space of the dynamical system for the = 
case. This figure is obtained by integrating the nonfinear sys- 
tem for initial conditions close to the fixed points. The (red) 
point on the left represents the unstable static solution while 
the (blue) point on the right represents the stable static so- 
lution. The continuous (red) curves intersecting the unstable 
point are parts of the separartix. 



two solutions collapse on the same point and disappear 
exploiting a saddle-center bifurcation. 



6 




(a) (b) 




H H 



(c) (d) 

FIG. 7: Dynamical behaviour of the system in the 
(H, p)— plane for initial conditions close to the static solutions. 
The continuous (red) line is a separatrix, a curve which crosses 
the unstable (upper, red) fixed point, enclosing the stable re- 
gion around the neutrally stable (lower, blue) fixed point. The 
phase space portraits are drawn using always the same arbi- 
trary values for the parameter and varying the the equation 
of state parameter, i.e. (a) w — —0.88, (b) w = —0.87, (c) 
w = -0.86667, (d) w = -0.865. 



V. CONCLUSIONS 

We have considered static solutions in the cosmologi- 
cal sector [21] of a recently proposed theory of massive 
gravity [TS] . We have shown that the effect of a massive 
graviton is to enrich the phase space of the cosmologi- 
cal equations, enlarging the ranges of existence of static 
solutions and affecting their stability properties. The so- 
lutions found exhibit stability properties rather different 
from those of the standard ES solution of GR, which 
requires a positive cosmological constant and a positive 
spatial curvature in addition to a perfect fluid with an 
equation of state parameter w > —1/3. 

Two kinds of solutions are present: neutrally stable so- 
lutions and unstable solutions (of the saddle type). Sum- 
ming up, for spatially closed (/C = 1) models three cases 
are possible: i) both the unstable and the neutrally sta- 
ble solutions are admitted, ii) either the unstable solu- 
tion or the neutrally stable solution is admitted, or in) no 
static solutions are admitted. For spatially flat (/C = 0) 
and open (/C = —1) models, three cases are possible: i) 
both the unstable and the neutrally stable solutions are 
admitted, ii) only the unstable solution is admitted, or 
Hi) no static solutions are admitted. Notice that, in the 



considered region of the parameter space, the neutrally 
stable solution requires a negative equation of state pa- 
rameter w; in particular it must be —1 < w < —2/3 in 
the /C = 0, — 1 cases and —1.2 < w < in the K, = 1 case. 

Our result implies the possibility of constructing mod- 
els in which the Universe oscillates indefinitely about an 
initial static state, thus the fine-tuning problem suffered 
by the emergent Universe scenario in GR [S] is amelio- 
rated when MG modifications are taken into account. On 
the other hand, this result raises the question of finding a 
mechanism in order to break the regime of infinite oscil- 
lations able to enter the current expanding phase under- 
gone by the Universe [7| . This result can be achieved by 
varying one of the model parameters, namely the equa- 
tion of state parameter w and the dimensionless param- 
eter C4 due to MG modification, in order for the system 
to undergo a bifurcation thus chancing the topological 
structure of the phase space. Such a mechanism has no 
dynamical explanation thus it looks quite unsatisfactory 
as well. Moreover, a full-fledged stability analysis against 
nonhomogeneous and/or nonisotropic modes would prob- 
ably reveal other instabilities as it happens in GR but 
such an analysis is well beyond the scope of this paper. 
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Lyapunov stable but not attracting [26] 
A static solution with open spatial sections is admitted 
in GR, it is unstable of the saddle type and requires very 
peculiar conditions, namely, a negative cosmological con- 
stant and a phantom fluid [14] 



